module setthmod
	use basicmodule
	implicit none

	integer				:: nthconst=4
	
	real				:: thq=0.90
	
 	real				:: e1ekr=1.2
	
	integer, parameter	:: ntailspecial=1, ntailbound=2
	integer				:: tailflag(2), tfcounter=0

	integer				:: ii
	integer, parameter	:: nxsi=25
	real, parameter		:: xsiminc=-1
	real, parameter		:: xsigrid(nxsi)=[(xsiminc+(xsimax-xsiminc)*(ii-1)/(nxsi-1),ii=1,nxsi)]
	real				:: muppergrid(nxsi),maxxxflat
	
	contains
		
	subroutine setmupper
		integer, parameter	:: nsim=20000
		real	:: xsi, eps(k,nsim), Z(k), zs(k,nsim),ml,dm,m,meas(nsim),b(2)
		integer	:: ix,il,l,i
		
		do l=1,nsim
			call rnun(eps(:,l))
		enddo
		do ix=1,nxsi
			xsi=xsigrid(ix)
			if(abs(xsi)<1E-6) xsi=1E-6
			do l=1,nsim
				z=-log(eps(:,l))
				do i=2,k
					z(i)=z(i-1)+z(i)
				enddo
				zs(:,l)=(z**(-xsi)-1)/xsi
			enddo
			ml=getmlower(xsi)
			dm=(expecval(1,xsi)-1.1*expecval(k,xsi))/(1.1-1)
			do i=1,20
				m=ml+dm
				if(count((m+zs(1,:))/(m+zs(k,:))<swflat)<thq*nsim) ml=m
				dm=0.5*dm
			enddo
			muppergrid(ix)=m
		enddo
		ml=0
		dm=1
		do i=1,20
			xsi=xsimin+(xsimax-xsimin)*(ml+dm)
			if(getmupper(xsi)<1/xsi .or. xsi<0) ml=ml+dm
			dm=0.5*dm
		enddo
		maxxxflat=ml
		print *,"largest xsi for flatness",xsimin+(xsimax-xsimin)*maxxxflat 
				
	end subroutine
	
	function getmupper(xsi) result(val)
		use qdval_int
		real	:: xsi,val
		integer	:: ialt
		val=qdval(xsi,xsigrid,muppergrid,check=.false.)
	end function

	function getsiglower(m,xsi,t) result(val)	 ! Y(1) takes on values smaller than t thq percent of the time
		real	:: xsi,m,val,t
		real	:: q
		q=(-1 + m*xsi + (-Log(thq))**(-xsi))/xsi
		val=t/q
	end function
	
	subroutine settailflag
		integer	:: i,j
		select case(stage)
			case(0)
				if(tfcounter<=ntailspecial**2) then
					tailflag(1)=-(mod(tfcounter-1,ntailspecial)+1)
					tailflag(2)=-((tfcounter-1)/ntailspecial+1)
				else
					tailflag(1)=mod(tfcounter,ntailbound)+1
					j=mod(tfcounter/ntailbound,10)+1
					if(j<=ntailspecial) then
						tailflag(2)=-j
					else
						tailflag(2)=mod(j,2)+1
					endif
				endif
			case(1)
				j=mod(tfcounter,10)+1
				if(j>ntailspecial) then
					tailflag=[0,mod(j,ntailbound)+1]
				else
					tailflag=[0,-j]
				endif
			case(2)
				tailflag=[0,0]
			case(3)
				tailflag(1)=mod(tfcounter,2)+10
				j=mod(tfcounter/2,10)+1
				if(j<=ntailspecial) then
					tailflag(2)=-j
				else
					tailflag(2)=mod(j,2)+10
				endif
			end select
			tfcounter=tfcounter+1
	end subroutine
	
	function getthfromx(x) result(val)
		real	:: x(6),val(6)
		val=[get_msx0t(x(1:3),tailflag(1)),get_msx0t(x(4:6),tailflag(2))]
	end function
	
	subroutine setth0(x,th)
		real	:: x(6),th(6)
		th=getthfromx(x)
	end subroutine
	
	function gets2sum(m,xsi,kbar) result(val)
		real	:: xsi,m,val,mu
		integer	:: i,kbar
		real	:: aj1,aj2
		val=0		  
		do i=k+1,kbar
			val=val+(Gamma(i-2*xsi)+2*(xsi*m-1)*Gamma(i-xsi))/Gamma(real(i))+(xsi*m-1)**2
		enddo
		val=val/xsi**2
		return
		if((expecval(kbar,xsi)+m)<0) then
			mu=kbar*(gamma(1+kbar-xsi)/((1-xsi)*gamma(kbar+1.0))+(m*xsi-1))/xsi
			if(mu<0) then
				print *,"mu<0"
			endif
			val=val+abs(expecval(kbar,xsi)+m)*mu
		endif
	end function
	
	function getsigupper(m,xsi) result(val)
		real	:: xsi,m,val
		integer	:: alt
		real	:: ekbarm
		val=gets2sum(m,xsi,n0-k)
		val=val**(-0.5)
		val=sqrt(varestfac)*val
	end function


	function get_msx0t(x0,tflag) result(msx)
		integer	:: tflag
		real	:: x0(3),x(3),msx(3)
		real	:: xsi, sig, m, b(2), s2
		real, parameter	:: e1ekr=1.05

		select case(tflag)
			case(0)
				msx=get_msx0(x0)
			case(-1)								! tiny tail				
				msx=99
			case(-2)								! large but flat
				xsi=xsimax
				m=(expecval(1,xsi)-e1ekr*expecval(k,xsi))/(e1ekr-1)
				sig=getsigupper(m,xsi)
				msx=[m,sig,xsi]
			case(-3)
				xsi=xsimin
				m=(expecval(1,xsi)-e1ekr*expecval(k,xsi))/(e1ekr-1)
				sig=getsigupper(m,xsi)
				msx=[m,sig,xsi]
			case(1)						! flat
				x=x0
				x(3)=x(3)*maxxxflat
				x(1)=1
				msx=get_msx0(x)
			case(2)						! small
				x=x0
				x(2)=0
				msx=get_msx0(x)
			case(3)						! small
				x=x0
				x(3)=1
				x(1)=1
				x(2)=x(2)-1
				msx=get_msx0(x)
			case(10)					! one tail for check
				x=x0
				msx=get_msx0(x)		
			case(11)
				x=x0
				msx=get_msx0(x)			! other tail for check

		end select
	end function
	
	function getmlower(xsi) result(val)
		real	:: xsi,val,w,mu
		val=-(gamma(1+n0-xsi)/((1-xsi)*gamma(n0+1.0))-1)/xsi		
	end function
	
	function get_msx0(x0) result(msx)
		real	:: x0(3),x(3),msx(3)
		real	:: xsi, sig, m, b(2), s2
		x=x0
		xsi=xsimin+(xsimax-xsimin)*x(3)
		if(abs(xsi)<1E-6) xsi=1E-6*merge(1,-1,xsi>0)
		b(1)=getmlower(xsi)
		b(2)=min(getmupper(xsi),merge(1/xsi,10000.0,xsi>0))
		b(2)=max(b(1),b(2))
		m=b(1)+(b(2)-b(1))*x(1)
		b(1)=getsiglower(m,xsi,swsmall)
		b(2)=getsigupper(m,xsi)
		b(2)=max(b(1),b(2))
		sig=b(1)*(b(2)/b(1))**x(2)
		
		msx=[m,sig,xsi]
	end function

	function getth_const(i,th) result(val)
		real	:: th(nth),val
		integer	:: i,j
		real	:: m
		select case(i)
			case(1)
				val=getej(n0/2,th(1:3))+getej(n0/2,th(4:6))
			case(2)
				m=0
				do j=1,n0/2
					m=m+getej(j,th(1:3))-getej(j,th(4:6))
				enddo
				val=max(m,0.0)*getej(n0/2,th(4:6))
			case(3)
				m=0
				do j=1,n0/2
					m=m+getej(j,th(4:6))-getej(j,th(1:3))
				enddo
				val=max(m,0.0)*getej(n0/2,th(1:3))
			case(4)
				val=varestfac-th(2)**2*gets2sum(th(1),th(3),n0/2)-th(5)**2*gets2sum(th(4),th(6),n0/2)
		end select
	end function
	
	function is_th_valid(th) result(val)
		real	:: th(nth)
		logical	:: val
		integer	:: i
		val=.true.
		if(any(th==99)) return
		do i=1,nthconst
			if(getth_const(i,th)<-1E-3)	then
				val=.false.
				return
			endif
		enddo
	end function
	
	
	function getf1(Y) result(val)
		real	:: Y(k), val
		real	:: xs(k),s,xsi
		integer	:: i,j
		integer, parameter	:: nxsigrid=10
		real, parameter		:: xsigrid(nxsigrid)=[(1E-2+xsimin1+(xsimax1-xsimin1)*(i-1.0)/(nxsigrid-1),i=1,nxsigrid)]
		real, parameter		:: q=1.0
		real	:: wm
		if(y(1)/y(k)<1.05 .and. y(k)>0) then
			val=1E100
			return
		endif
		xs=(y-y(k))/(y(1)-y(k))
		val=0
		do j=1,nxsigrid
			xsi=xsigrid(j)
			wm=1 !exp(xsi)
			do i=1,nGQ
				s=(GQxw(i,1)**(-xsi)-1)/xsi
				val=val+wm*GQxw(i,2)*gamma(k-q*xsi)*GQxw(i,1)**(-xsi-1)*exp(-(1+1/xsi)*sum(log(1+xsi*s*xs))+(k+q-2)*log(s))
			enddo
		enddo
		val=val/nxsigrid
		val=val/(y(1)-y(k))**(k+q-1)	
	end function

end module

